High-Precision c and b Masses, and QCD Coupling 
from Current-Current Correlators in Lattice and Continuum QCD 
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We extend our earlier lattice-QCD analysis of heavy-quark correlators to smaller lattice spac- 
ings and larger masses to obtain new values for the c mass and QCD coupling, and, for the first 
time, values for the b mass: mc(3GeV,n/ = 4) = 0.986(6) GeV, a'^(Mz,nf = 5) = 0.1183(7), 
and 7715(10 GeV, n/ = 5) = 3.617(25) GeV. These are among the most accurate determinations 
by any method. We check our results using a nonperturbative determination of the mass ra- 
tio mb{i-i,nf)/mc{j-L,nf); the two methods agree to within our 1% errors and taken together im- 
ply rrib/mc =4.51(4). We also update our previous analysis of q,^ from Wilson loops to account 
for revised values for ri and ri/a, finding a new value Q]^(Mz, n/ = 5) = 0.1184(6); and we update 
our recent values for light-quark masses from the ratio mc/nis. Finally, in the Appendix, we derive 
a procedure for simplifying and accelerating complicated least-squares fits. 

PACS numbers: 11.15.Ha,12.38.Aw,12.38.Gc 
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I. INTRODUCTION 

Precise values for the QCD coupling ajj^ and the quark 
masses are important for high-precision tests of the Stan- 
dard Model of particle physics. In a recent paper we 
showed how to use realistic lattice QCD simulations to 
extract both the coupling and the charm quark's mass TOc 
from zero-momentum moments of correlators built from 
the c quark's (UV cutoff- independent) pseudoscalar den- 
sity operator mc'0c75'0c [!]• In this paper we refine our 
previous analysis and extend it to include other quark 
masses, up to and including the 6-quark mass. As a re- 
sult our coupling constant and mass determinations from 
these correlators are among the most accurate by any 
method. 

Low moments of heavy-quark correlators are perturba- 
tive and several are now known through Ofa^) in pertur- 
bation theory (that is, four-loop order) Moments of 
correlators built from the electromagnetic currents can be 
estimated nonperturbatively, using dispersion relations, 
from experimental data for the electron-positron annihi- 
lation cross section, ^(e+e" ^-7* ^-X). Accurate values 
for both the c and b masses can be obtained by comparing 
these perturbative and nonperturbative determinations 
of the moments (for a recent discussion see [7]). 

In our earlier paper we showed that heavy-quark corre- 
lator moments are easily and accurately computed non- 
perturbatively using lattice QCD simulations, in place 
of experimental data, provided; 1) the electromagnetic 
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current is replaced by the pseudoscalar density multi- 
plied by the bare quark mass; 2) the discretization of the 
quark action has a partially conserved axial vector cur- 
rent (PCAC); and 3) the discretization remains accurate 
when applied to heavy quarks. In our simulations we 
use the HISQ discretization of the quark action, which 
is a highly corrected version of the standard staggered- 
quark action It has a chiral symmetry (PCAC) and 
has been used in a wide variety of accurate simulations 
involving c quarks [sl-fl^. 

Here we show that the HISQ action can be pushed to 
still higher masses — indeed, very close to the b mass — 
on new lattices, from the MILC collaboration [l^, with 
the smallest lattice spacings available today (0.044fm). 
Currently most high-precision lattice work on b physics 
relies upon nonrelativistic effective field theories, like 
NRQCD [iS El [il El. In this paper we show how 
to obtain accurate b physics using the fully relativistic 
HISQ action on these new lattices. 

In what follows, we first review how the QCD coupling 
and quark masses are extracted from heavy-quark cor- 
relators, in Section II. Then in Section HI we describe 
our lattice QCD simulations and discuss in detail the 
chief systematic errors in our simulation results. In Sec- 
tion IV we describe our fitting procedure and the re- 
sults of our analysis of the heavy-quark correlators. We 
check our calculation using a different, nonperturbative 
method to determine mb/rric in Section V. We then, in 
Section VI, update our previous calculation of the QCD 
coupling from Wilson loops to compare with our new re- 
sult from the correlators. We summarize our findings in 
Section VII and compare our results with work by oth- 
ers. There we also update our recent calculations of the 
light-quark masses from the c mass. In the Appendix we 
present a powerful simplification for complicated least- 
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squares fits that can greatly reduce the computing re- 
quired for fits. We use this technique in dealing with 
finite-a errors in our analysis. 



II. HEAVY-QUARK CORRELATOR MOMENTS 

Following our earlier paper , we focus on correlators 
formed from the pseudoscalar density of a heavy quark, 

j5=i'hl5ll'h-- 



Git) = a'Y^iamohf (0| j5(x, t)j,iO, 0)|0) 



(1) 



where moh is the heavy quark's bare mass (from the lat- 
tice QCD lagrangian), t is euclidean and periodic with 
period T, and the sum over spatial positions x sets the 
total three momentum to zero. In our earlier paper we 
examined only c quarks; here we will consider a range 
of masses between the c and b masses. While we have 
written this formula for use with the lattice regulator, 
it is important to note that the correlator is UV-finite 
because we include the factors of amoh- Consequently 
lattice and continuum G{t)s are equal when i 7^ up to 
©((am^)™) corrections, which vanish in the continuum 
limit. 

The moments of G(t) are particularly simple to ana- 
lyze: 



G„ = ^(i/«rG(i) 



(2) 



where, on our periodic lattice, 

t/ae {0, 1, 2... T/2a - 1, 0, -r/2a +1...-2, -1}. (3) 

Low moments emphasize small ts and so are perturbative; 
and moments with 71 > 4 are UV-cutoff independent. 
Therefore 



Gn - — 7 . ^^„_4 1- G>{{amh 



(4) 



for small n > 4, where m/i(/i) is the heavy quark's MS 
mass at scale /i, and the dimensionless factor gn can be 
computed using continuum perturbation theory. 

Again following our previous paper, we introduce re- 
duced moments to suppress both lattice artifacts and 
tuning errors in the heavy quark's mass 16]: 



R„ 



'"iih 



2amoh 



l/(n-4) 



for n — A, 



for 71 > 6, 



(5) 



where gI"'' is the moment in lowest-order, weak-coupling 
perturbation theory, using the lattice regulator, and m,,^ 
is the (nonperturbative) mass of the pseudo-Goldstone 
hh boson. The reduced moments can again be written in 
terms of continuum quantities: 



Rr 



for n 



z{fi/mh,mjjjrn{a^,fj,/mh) for n > 6, 



(6) 



up to 0{{amji)'^as) corrections, where 



2m/i(/z)' 



(7) 



and r„ is obtained from gn (Eq. (g])) and its value, gn\ 
in lowest-order continuum perturbation theory: 



l/(n-4) 



for n = 4, 
for n > 6. 



(8) 



Our strategy for extracting quark masses and the QCD 
coupling relies upon lattice simulations to determine non- 
perturbative values for the using simulation results 
for am,^^/ aniQh- We then compare this simulation "data" 
to the continuum perturbation theory formulas (Eq. ([6])). 
That is, we find values for aj^(/i) and z{^/mh, rari^) that 
make lattice and continuum results agree for small n > A. 
The function z{^/mhTmrif^) can then be combined with 



experimental results for tti^^ and rriri^ 
for the c and b quarks: 



to obtain masses 



2z{^./mc,7n^J') 



nibifJ.) 



Vb 



2z{^i/mb,nh^^) 
(9) 

Parameter ji sets the scale for in the perturbative 
expansions of the r„ . An obvious choice for this param- 
eter is /i = vrih since the quark mass, together with n, 
sets the momentum scale in our correlators. As noted 
in our previous paper, however, perturbation theory is 
somewhat more convergent if we use larger /is in the c- 
quark case. Consequently here we take /i/TOft = 3, which 
is approximately what we did in our previous paper. 

The mass and coupling determinations were done sep- 
arately in our previous paper. Here we extract them si- 
multaneously, to guarantee consistency between results. 
Also in our previous paper we considered only heavy- 
quark masses near the c mass. Here we explore a variety 
of masses ranging from just below the c mass to just 
below the b mass. This allows us to obtain a value for 
6-quark's mass. 



III. LATTICE QCD SIMULATIONS 
A. Simulation Results 

The gluon-configuration sets we use were created by 
the MILC collaboration. The relevant simulation param- 
eters are listed in Table HI 

Given a lattice spacing, the QCD action is specified 
completely by the values of the bare coupling constant 
and the bare quark masses. In our analyses we set the 
u and d quark masses equal; this approximation results in 
negligible errors (<C 1%) for the quantities studied in this 
paper. It is too costly to simulate QCD at the correct 
value for the u/d mass; we typically use masses that are 
2-5 times too large and extrapolate to values that give 
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TABLE I: Parameter sets used to generate the gluon con- 
figurations analyzed in this paper. The lattice spacing is 
specified in terms of the static-quark potential parameter 
ri = 0.3133(23) fm; values for ri/a are from [l^. The bare 
quark masses are for the ASQTAD formalism and uq is the 
fourth root of the plaquette. The spatial (L) and temporal 
(T) lengths of the lattices are also listed, as are the number 
of gluon configurations (Nd) and the number of time sources 
(A'ts) per configuration used in each case. Sets with similar 
lattice spacings are grouped. 



Set 


ri/a 


auorriou/d 




no 


L/a 


T/a 


iVcf X TVts 


1 


2.152(5) 


0.0097 


0.0484 


0.860 


16 


48 


631 X 2 


2 


2.138(4) 


0.0194 


0.0484 


0.861 


16 


48 


631 X 2 


3 


2.647(3) 


0.005 


0.05 


0.868 


24 


64 


678 X 2 


4 


2.618(3) 


0.01 


0.05 


0.868 


20 


64 


595 X 2 


5 


2.618(3) 


0.01 


0.05 


0.868 


28 


64 


269 X 2 


6 


3.699(3) 


0.0062 


0.031 


0.878 


28 


96 


566 X 4 


7 


3.712(4) 


0.0124 


0.031 


0.879 


28 


96 


265 X 4 


8 


5.296(7) 


0.0036 


0.018 


0.888 


48 


144 


201 X 2 


9 


7.115(20) 


0.0028 


0.014 


0.895 


64 


192 


208 X 2 



the correct mass for the Tr^-meson. We tune the strange 
quark mass to give the correct mass for the (fictitious) 
rys meson [lo| . The c and b masses are tuned to give 
correct masses for the 77c and 77b mesons, respectively. 

It is convenient in QCD simulations to specify a value 
for the bare coupling constant and then extract the value 
of the lattice spacing from the simulation. We set the 
lattice spacing using MILC results for ri/a, computed 
from the heavy-quark potential, and [loj 

n = 0.3133(23) fm. (10) 

The MILC configurations include vacuum polarization 
contributions from only the lightest three quark flavors, 
using the ASQTAD discretization. Vacuum polarization 
effects from the heavier c and b quarks are easily incor- 
porated into our final results for quark masses and the 
QCD coupling using perturbation theory. 

We computed heavy-quark correlators (Eq. (jlj) using 
the HISQ discretization p] for a variety of bare heavy- 
quark masses amoh on the MILC gluon configurations. 
Our results for the reduced moments i?„ with n = 4-18 
are given in Table |IT1 

In Table |lTl we also give masses am^^ from the sim- 
ulations for the pseudo-Goldstone meson made from 
two heavy quarks. These were computed using single- 
exponential fits to G{t) for the middle 30% of ts on the 
lattice for all configurations except the two smallest lat- 
tice spacings where we used only 8% of the ta. We have 
less statistics for the two finest lattice spacings and con- 
sequently the fits did not work as well for these. We 
increased the statistical errors on our results for am^^ by 
factors of 1.4 and 2 for the next-to- finest and finest lat- 
tice spacings (sets 8 and 9), respectively, to account for 
this. The statistical errors here are very small and have 
only a small impact on our final results. We also verified 
our results with multi-exponential fits in every case. 



B. Systematic Errors 

As discussed above, our goal is to find values for 
a^(/i) and z{iJi/mh,ra^^) (Eq. (O) that make the the- 
oretical results from perturbative QCD agree, to within 
statistical and systematic errors, with Monte Carlo sim- 
ulation "data" for the reduced moments. We simultane- 
ously analyze results for all of our lattice spacings and 
most of our masses, and for moments with 4 < rt < 10. 
We focus on these particular moments for our final re- 
sults since their perturbation theory is known to third 
order. 

Systematic errors are larger here than statistical errors, 
which contribute less than 0.3%. We discuss the most 
important sources of systematic error in this section. 

1. ruh Extrapolations 

We need the m,,^ dependence of the mass-ratio func- 
tion z{fj,/mh = -i,'mri^) in order to extract c and b masses 
from our simulation (using Eq. (|9])). We parameterize 
this dependence as follows: 

^' / 2A V 
z(/i/m,,,m^J = ^Zj(^/TOft) , (11) 

where the ZjS are determined in our fit. This is an ex- 
pansion in the QCD scale, which we take to be 

A = 0.5GeV, (12) 

divided by to^^/2, which we use as a proxy for the quark 
mass. The expansion is adequate for the range of quark 
masses used in our analysis, where (2A/m^^)^ ranges ap- 
proximately between l/rn^^ = 0.01 and (1/m^J^ — 0.1; 
the singular point rrih = is infinitely far away in this 
parameterization. In our fits we keep terms only through 
order iV^ = 4, but, as we discuss later, our results are un- 
changed by additional terms. On dimensional grounds, 
we assume a priori that the coefficients are 

z,(3) = 0±l. (13) 

2. Finite- Lattice Spacing Errors 

Discretization errors are of order {arrihf'^as for i>l. 
We model these by 

= i?„(/i, m^,, a, Na^), (14) 

where: fit function Rn(p,mj^^,a,Nam) has the double 
expansion 

Rn{ti,m^^,a,Nam) = RT'/ (15) 
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TABLE II: Results for the reduced moments Rn and pseudoscalar-meson mass am^^ obtained from (n/ = 3) simulations using 
different bare heavy-quark (HISQ) masses amoh and gluon conhguration sets (see Table [I]). The errors listed here are statistical 
errors from the Monte Carlo simulation. Results where am^^ > 1.95 are omitted from our final analysis, as are R„s with n> 10. 



Set 


amoh 




i?4 


Re 


Rs 




Rio 




Rl2 




Rl4 




Rie 


Ris 




1 


0.660 


1 riono 


1 \ 
i j 


1.2132(3) 


1.5364(3) 


1.4151 


2) 


1.3476 


1) 


1.3001 


1) 


1.2649 


1) 


1.2378(1) 


1.2164 


1) 




0.810 


z.iyoo 


i ) 


1.1643(2) 


1.4427(2) 


1.3619 


1) 


1.3148 


1) 


1.2780 


1) 


1.2481 


1) 


1.2238(1) 


1.2039 


1) 




0.825 




i ) 


1.1604(2) 


1.4339(2) 


1.3563 


1) 


1.3111 


1) 


1.2754 


1) 


1.2462 


1) 


1.2222(1) 


1.2025 


1) 


2 


0.825 


z.ziyt) 


i ) 


1.1591(2) 


1.4327(2) 


1.3556 


1) 


1.3106 


1) 


1.2751 


1) 


1.2459 


1) 


1.2221(1) 


1.2024 


1) 


3 


0.650 


1 Q/1 CO 


1 j 


1.1809(2) 


1.4805(2) 


1.3755 


1) 


1.3160 


1) 


1.2740 


1) 


1.2429 


1) 


1.2190(1) 


1.2000 


1) 


4 


0.440 




1 \ 
i ) 


1.2752(4) 


1.6144(4) 


1.4397 


2) 


1.3561 


2) 


1.3041 


1) 


1.2678 


1) 


1.2408(1) 


1.2200 


1) 




0.630 


1 or»o c 


1 ^ 
i) 


1.1881(3) 


1.4935(2) 


1.3826 


1) 


1.3205 


1) 


1.2773 


1) 


1.2456 


1) 


1.2214(1) 


1.2021 


1) 




0.660 


1 Qfifi'7 


i ) 


1.1782(2) 


1.4764(2) 


1.3738 


1) 


1.3152 


1) 


1.2736 


1) 


1.2426 


1) 


1.2187(1) 


1.1997 


1) 




0.720 


i.yoii 


i ) 


1.1605(2) 


1.4435(2) 


1.3559 


1) 


1.3044 


1) 


1.2662 


1) 


1.2367 


1) 


1.2136(1) 


1.1950 


1) 




0.850 


01 0/1 


i j 


1.1301(2) 


1.3763(1) 


1.3145 


1) 


1.2774 


1) 


1.2473 


1) 


1.2221 


1) 


1.2012(1) 


1.1839 


1) 


5 


0.630 


1 Qf\Qfi 
i.oUOD 


i ) 


1.1882(1) 


1.4936(1) 


1.3826 


1) 


1.3205 


1) 


1.2774 


1) 


1.2457 


1) 


1.2214(0) 


1.2022 


0) 


6 


0.300 


i.Uol4 


1 \ 
i ) 


1.2930(3) 


1.6061(3) 


1.4249 


2) 


1.3444 


1) 


1.2953 


1) 


1.2610 


1) 


1.2353(1) 


1.2153 


1) 




0.413 


i.zoUD 


i ) 


1.2224(2) 


1.5216(2) 


1.3796 


1) 


1.3115 


1) 


1.2689 


1) 


1.2390 


1) 


1.2164(1) 


1.1985 


1) 




0.430 


i.olDy 


i ) 


1.2145(2) 


1.5113(2) 


1.3743 


1) 


1.3076 


1) 


1.2658 


1) 


1.2363 


1) 


1.2141(1) 


1.1964 


1) 




0.440 


1 OQOO 


ij 


1.2100(2) 


1.5054(2) 


1.3712 


1) 


1.3054 


1) 


1.2640 


1) 


1.2348 


1) 


1.2127(1) 


1.1952 


1) 




0.450 


i.ooyo 


i ) 


1.2057(2) 


1.4996(2) 


1.3683 


1) 


1.3033 


1) 


1.2623 


1) 


1.2333 


1) 


1.2114(1) 


1.1941 


1) 




0.700 


1 QfiKA 


1 ^ 
1 } 


1.1301(1) 


1.3782(1) 


1.3053 


1) 


1.2616 


1) 


1.2294 


1) 


1.2048 


0) 


1.1857(0) 


1.1705 


0) 




0.850 


z.i4yo 


i ) 


1.1026(1) 


1.3163(1) 


1.2671 


1) 


1.2366 


0) 


1.2114 


0) 


1.1903 


0) 


1.1729(0) 


1.1584 


0) 


7 


0.427 


i.oU / 4 


i ) 


1.2131(3) 


1.5091(3) 


1.3729 


2) 


1.3066 


1) 


1.2651 


1) 


1.2358 


1) 


1.2137(1) 


1.1961 


1) 


8 


0.273 


0.8993 


3) 


1.2454(8) 


1.5234(9) 


1.3739 


7) 


1.3069 


6) 


1.2657 


6) 


1.2366 


6) 


1.2145(6) 


1.1969 


5) 




0.280 


0.9154 


2) 


1.2403(5) 


1.5175(5) 


1.3706 


3) 


1.3045 


3) 


1.2638 


2) 


1.2350 


2) 


1.2132(2) 


1.1958 


2) 




0.064 


1.5254 


1) 


t.i324(2j 


i.36(4(2j 


1 OO 

i.200/ 


2) 


i.240o 


i j 


1 oi no 
i.ziOz 


i j 


i.isoo 


t j 


t.i / i9(ij 


1 1 

i.ioo< 


i j 




0.705 


1.8084 


1) 


1.1043(2) 


1.3156(2) 


1.2574 


1) 


1.2217 


1) 


1.1952 


1) 


1.1750 


1) 


1.1593(1) 


1.1467 


1) 




0.760 


1.9157 


1) 


1.0955(2) 


1.2965(2) 


1.2460 


1) 


1.2142 


1) 


1.1895 


1) 


1.1701 


1) 


1.1547(1) 


1.1423 


1) 




0.850 


2.0875 


1) 


1.0831(2) 


1.2666(1) 


1.2266 


1) 


1.2010 


1) 


1.1799 


1) 


1.1621 


1) 


1.1474(1) 


1.1353 


1) 


9 


0.195 


0.6710 


2) 


1.2583(5) 


1.5243(5) 


1.3733 


4) 


1.3066 


3) 


1.2655 


3) 


1.2364 


2) 


1.2144(2) 


1.1968 


2) 




0.400 


1.1325 


2) 


1.1532(3) 


1.3800(3) 


1.2838 


2) 


1.2370 


2) 


1.2077 


2) 


1.1869 


2) 


1.1710(2) 


1.1583 


2) 




0.500 


1.3446 


2) 


1.1267(2) 


1.3410(2) 


1.2616 


1) 


1.2198 


1) 


1.1927 


1) 


1.1734 


1) 


1.1588(1) 


1.1471 


1) 




0.700 


1.7518 


1) 


1.0900(1) 


1.2765(1) 


1.2261 


1) 


1.1949 


1) 


1.1718 


1) 


1.1542 


1) 


1.1407(1) 


1.1299 


1) 




0.850 


2.0428 


1) 


1.0712(1) 


1.2327(1) 


1.1983 


1) 


1.1760 


1) 


1.1574 


1) 


1.1418 


1) 


1.1290(1) 


1.1185 


1) 



the c^J'^s are determined in our fit, is given 



by Eq. 



i+j< niax(Afam,^z), 



(16) 



and again we use m,j,^/2 in place of the quark mass. 
This expansion allows for finite-a corrections involving 
(am,,^/2)^, (aA)^, and cross terms, with m,,^ -dependent 
coefficients. We assume a priori that 



c|;' = o±2/n 



(17) 



which implies smaller a dependence for larger ns. This is 
expected (and obvious in our simulation data) since the 
reduced moments become more infrared as n increases. 
The exact functional form of the n dependence has little 
effect on our results, as we show later. 

In our fits we take Nz = 4. While low orders suffice 
for the 2A/m^^ expansion, expansion parameter aTO^^/2 
ranges between 0.3 and 1.1, and higher orders are nec- 
essary, especially given our tiny statistical errors. We 
find that our fit results don't converge well unless Nam is 
larger than 10-20. Also we have difficulty getting good 
fits if we include data with am^^ > 1.95 from Table HIl 



The amri,^/2 expansion may not converge for these last 
cases and therefore we exclude such data from our final 
analysis. 

The fit function has many more fit parameters c'"' 
than we have simulation data points when Nam is so 
large. This does not cause problems in (Bayesian) con- 
strained fits since the parameters' priors (Eq. (IT71) ) are 
included in the fit as extra data Each parameter 

has a prior and therefore we always have more data than 
parameters. 

It is, however, very time consuming to fit a function 
with so many fit parameters. Although it is not essential 
for our analysis, there is a trick that greatly accelerates 
this kind of fit. The idea is to fit a modified moment 
^Jf** in place of i?Jf" where 



latt 



R 



olatt I 
+ 

Nam 

latt 



(18) 



2A 



The modified moment is fit with the 
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much simpler formula (simpler since Nam Nam) 

R':'' ^Rn{ti,m^„a,Nam)- (19) 

where _R„(...) is again given by Eq. (|15p . To evalu- 
ate ^Jf" from Eq. ([T^ . we treat the coefficients c[j^ 
with i > Nam SiS new data with means and standard de- 
viations specified by the prior, Eq. (|17p . Uncertainties 

coming from the c^"^s are combined in quadrature with 
the statistical error in i?Jf** to obtain a new error esti- 
mate for ^Jf (but leaving the central value unchanged) . 
In effect we are increasing the error in the reduced mo- 
ment to account for high-order (am^^/2)^* terms omitted 
from the fit formula Eq. p^ . By choosing Nam ^ Nam, 
most of the terms are incorporated into .Rjf** 

(Eq. (|18p ). where they are inexpensive, and relatively 
few end up in the fit function Rn{- ■ ■) (Eq. (|15|) ). where 
they add parameters to the fit and increase its cost. Note 
that the new errors introduce correlations between ^Jf**s 
computed with different lattice spacings or quark masses, 
since the same c^J^'s are used for all as and »7i,,^s. These 
correlations are important and need to be preserved in 
the fit. 

Our procedure, whereby terms are moved out of the 
fitting function and incorporated into new (correlated) 
errors in the Monte Carlo fit data, is generally useful. 
Somewhat remarkably, final fit results are completely (or 
almost completely) independent of the number of terms 
that are transferred when fits are linear (or almost lin- 
ear) in the associated parameters. (The general theorem 
from which this result follows is proven in the Appendix.) 
Consequently, in our analysis here, we can take Nam very 
large — say. Nam = 80 — and still have very fast fits by 
keeping Nam very small. With Nam = 80 we find, for ex- 
ample, that setting Nam = in _R„(. . .) (no terms) gives 
essentially identical results for our quark masses and cou- 
pling as setting TVam = 30 (140 terms), even though the 
latter fit requires 22 times more computing. We used 
this procedure, with Nam = 0, for most of our testing and 
development in this project. 

3. Truncated Perturbation Theory 
The perturbative part, 

Afpth 

''"(o^Ms, fJ-fmh) = 1+51 ^nJi^^/^h)aL-{^l), (20) 

of the reduced moments is known at best through third 
order. We present coefficients r„j through j = 3 in Ta- 
bleHIlll^; the values for ri = 4-10 are exact, while r„3 
is estimated for the others. In our fits we include higher- 
order terms by treating the coefficients of these terms as 
fit parameters with prior 

r„j(l)=0±0.5 (21) 



TABLE III: Perturbation theory coefRcients (n/ = 3) for 
Tn Q-Q. Coefficients are defined by r„ = 1-1- ^rijaL—{fi) 
for /i = mh(/x). The third-order coefficients are exact for 
4 < 71 < 10. The other coefRcients are based upon estimates; 
we assign conservative errors to these. 



n 


rnl 


r„2 


J"n3 


4 


0.7427 


-0.0577 


0.0591 


6 


0.6160 


0.4767 


-0.0527 


8 


0.3164 


0.3446 


0.0634 


10 


0.1861 


0.2696 


0.1238 


12 


0.1081 


0.2130 


0.1(3) 


14 


0.0544 


0.1674 


0.1(3) 


16 


0.0146 


0.1293 


0.1(3) 


18 


-0.0165 


0.0965 


0.1(3) 



for any coefficient that hasn't been computed in pertur- 
bation theory. We set A'^pth = 6 since then contributions 
from still higher orders should be less than 0.1% (and 
setting iVpth = 8 doesn't change our results) . 

The perturbative coefRcients for fj./mh = l (Table UlI)) 
are small and relatively uncorrelated from order-to-order. 
This is less true for fi/mh = 3, which is where we wish 
to work (see Section II), because of log(^/m/i)'" terms. 
In order to capture these effects, we use renormaliza- 
tion group equations to express the rnj{-i) coefficients 
(for all j < A^pth) in terms of the r„j(l) coefficients 
and \og{^ / rrifi) , and substitute the results from Table Hill 
for j < 3 and from the prior (Eq. ((2T|) ') for j > 3. 
This procedure generates (correlated) priors for the un- 
known coefficients at fi/mh — -i that properly account for 
renormalization-group logarithms. 



4- Q^MS" Evolution 

As discussed above, we fix the ratio of ii/mh{fJ-) in our 
analysis. This means that the renormalization scale fi 
varies over a wide range of values for the different m^s 
we use. The coupling constant aj^{fi) used in the per- 
turbative expansions for the r„s is specified at /i = 5 GeV 
by fit parameter ap, with prior 

ao = 0.20 ±0.01, (22) 

where 

ao = aMs(5 GeV, n/ = 3). (23) 

The prior corresponds to ayf^{Mz) =0.118(3) — a very 
broad range, which means that the prior has little im- 
pact on our final fit results. The coupling value at any 
scale /i 5 GeV is obtained by integrating (numerically) 
the QCD evolution equation for ctjjg^fj.) starting with 
value ao at scale 5 GeV. We use the MS beta function 
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through sixth order in Oj^jg, 
,2 ^Q'ms"(a^) 



= -/^o"Ms-/^i"Ms-/32«tis (24) 



MS 



6 

MS' 



where /3o . . . ^3 are known from perturbation theory and 
/34 is taken as a fit parameter with prior 



/34 = ± (T^ 



(25) 



where cr^ is the root-mean-square average oi ■ ■ ■ [la, 
[Tgj . We include this last term to estimate the uncertain- 
ties in our final results caused by unknown terms in the 
beta function. 

Our simulations include vacuum polarization effects 
from only the three lightest quarks. We use perturbation 
theory, together with the c and b masses that come out 
of our analysis, to incorporate vacuum polarization ef- 
fects from the heavier quarks into our final results for the 
masses and QCD coupling (using formulas from [g^, [2l| 
to add the c and b quarks at scales ii = mc and mb, re- 
spectively). 



5. Nonperturbative Condensates 

As discussed in our previous paper, nonperturbative ef- 
fects dominate the reduced moments when n is large. The 
dominant nonperturbative contribution, which is from 
the gluon condensate, is quite small, however, for the 
range of ns and quark masses we use here. We correct 
for it by replacing 



R 



latt 



R: 



latt 



1 



(26) 



where dn is computed to leading order in perturbation 
theory [23 with rrih = mh{mh,), which we approximate 
by m,,^/2.27. We take 



{asC^/Ti) = 0±0.012GcV^ 



(27) 



which covers the range of most current estimates [23j. 
The correction factor in Eq. (|26l) adds (slightly) to the 
error in i?Jf'* (and introduces new correlations between 
different moments, since the same {asC^ /tt) is assumed 
for every moment, lattice spacing and quark mass). 



6. Finite Volume Errors 

We expect small errors due to the fact that our simu- 
lation lattices are only about 2.5 fm across. We allow for 
the possibility of finite-volume errors by replacing 



olatt , r? 
it_ — > K, 



latt 



1 + /, 



th 



R 



,pth 



(28) 



where AR^^^ is the finite volume error in leading-order 
perturbation theory and 



/„ = ± 0.5. 



(29) 



The true finite-volume errors are expected to be 
smaller, because of quark confinement, than the pertur- 
bative errors that we use to model them here. We verified 
this by running two sets of simulations that were iden- 
tical except for the spatial volume (gluon configuration 
sets 4 and 5 in Table HI. The differences between the two 
simulations are smaller than our statistical errors, but 
the statistical errors are much smaller than our estimate 
above. Our error estimate here is very conservative, but 
has negligible impact on our final results. 



7. Sea-Quark Masses 

The sea-quark masses used in our simulations are not 
exactly correct. To correct for this we replace 



l + > 



2Smi + 6ms 



(30) 



where Smi and Snig are the errors in the u/d and s masses 
(see [13 for more details), respectively, and 



g„ = 0±0.01. 



(31) 



This correction introduces (correlated) errors into 
the i?Jf"s that are of order 0.5-1%. Direct comparison 
of results from configuration sets 6 and 7 (or 1-2 and 3- 
4) in Table H] suggests that sea-quark mass effects are no 
larger than 0.1%, so our error estimate is conservative. 

We have only included the leading dependence on the 
sea-quark mass, which comes from nonperturbative (chi- 
ral) effects. Quadratic terms from perturbation theory 
and other nonperturbative sources are negligible. 



IV. ANALYSIS AND RESULTS 

We have computed reduced moments for 30 different 
sets of lattice spacing, lattice volume and quark masses 
(Table |TT| . To extract quark masses and the QCD cou- 
pling, we fit moments with 4 < rt < 10 from 22 of these 
parameter sets (the ones with am,,^ < 1.95) — 88 pieces 
of simulation data in all. In this section we first describe 
the fitting method used to extract the masses and cou- 
pling, and then we review our results. 



A. Constrained Fits 

We analyze all four i?„s for all 22 parameter sets si- 
multaneously using a constrained fitting procedure based 



7 



upon Bayesian ideas [17| . In this procedure we minimize 
an augmented function of the form 



X 



where: 



(32) 



R 



latt 



the i?Jf** come from Table |TI] with corrections from 
Eqs. dM]), (EH) and dSO]); fit function i?„(. ..) is defined 
by Eq. (jisp: and tr^ is the error covariance matrix for 
the i?Jf". The sums i,j are over the 22 sets of lattice 
spacings and quark masses; the sums n, m range over of 
the moments 4, 6, 8, 10. 

Function R„(iJ.i,mri^i,ai, Nam) depends upon a large 
number of parameters, all of which are varied in the fit 
to minimize x^- Priors Sx^ are included for each of these: 



''Vh 



• parameters Zj, with prior Eq. (|13p . from the l/r 
expansion of z(fj,/mh, rriri^); 

• parameters c^"\ with prior Eq. (|17p . from the 
finite-lattice spacing corrections; 

• unknown perturbative coefficients with prior 
Eq. dUl) (evolved to ^/m^^B); 

• coupling parameter log(Q!o), with prior Eq. ([22]); 

• Pi in the QCD /3- function, with prior Eq. (^5]) : 

• lattice spacings for each giuon configuration set, 
with priors specified by simulation results for ri/a 
(Table H]) and the current value for ri (Eq. ([TU)) ): 

• values for arn^^i, with priors specified by our sim- 
ulation results (Table HI]) . 

The renormalization scales /i^ are obtained from the ratio 
fj./mh = 3, simulation results for to,,^, and Eq. ([7]). We 
take TVam = 30 for our final results. 



B. Results 



We fit our simulation data for the reduced mo- 
ments i?Jf" (Table HH) using fit function i?„(. . .) 



(Eq. (ITSt ) with Nam =30, as discussed in the previous 
section. The best-fit values for parameters zj give us the 
mass-ratio function z{^/mh = 3,rn,,^) (Eq. ([7])), which 
we plot in Figure [T] We also show our simulation re- 
sults there for i?^**/r„, together with the best- fit lines 
for each lattice spacing. Results are shown for the three 
moments that depend upon z, 5 different lattice spacings, 
and quark masses ranging from below the c mass almost 
to the b mass. The simulation data were all fit simulta- 
neously, using the same functions z(3,m^^) and aj^{^) 
(with /i = 3m^^/(2z)) for all moments. The fits are ex- 
cellent, with x^/88 = 0.19 for the 88 pieces of simulation 
data we fit. 
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FIG. 1: Function z(/^/mh = 3, m,,^) = m^^/(2mh) as a func- 
tion of m,,^ . The solid line, plus gray error envelope, shows 
the a = extrapolation obtained from our fit. This is com- 
pared with simulation results for R„/rn for n = 6, 8, 10 from 
our 5 different lattice spacings, together with the best fits 
(dashed lines) corresponding to those lattice spacings. Dashed 
lines for smaller lattice spacings extend further to the right. 
The points marked by an "x" are for the largest mass we 
tried (last line in Table |Il]); these are not included in the fit 
because am,,^ is too large. Finite-a errors become very small 
for the larger-n moments, causing points from diflterent lattice 
spacings to overlap. 



Evaluated at m^^ = 2.985(3) GeV |23|, the mass- 
ratio function is 2:(3,to^^) — 1.507(7). Combining this 
with Eq. ^ and perturbation theory, we can obtain the 
following results for the MS c-quark mass at different 
scales: 



TOc(3mc, n/ = 3) = 0.991(5) GeV, 
TOc(3GeV,n/ = 4) = 0.986(6) GeV, 
mc{mc, n/ = 4) = 1.273(6) GeV. 



(34) 



Similarly at m^^ = 9.395(5) GeV [23, the mass-ratio 
function is z(3, mrf^) = 1.296(8), and we obtain the follow- 
ing results for the MS 6-quark mass at different scales: 

m,,(3m,„ = 3) = 3.622(22) GeV. (35) 
mf,(10GeV,7i/ = 5) = 3.617(25) GeV, 
mb{mb,nf = 5) = 4.164(23) GeV. 

Note that the ratio mi,{ij,,nf)/mi.{^i,nf) is indepen- 
dent of fi and n f . We obtain the following result for this 
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FIG. 2: QCD coupling aj^(/i,n/ =3) as a function of m^^ 
where /i = 3mh. The solid line, plus gray error envelope, shows 
the best-fit coupling from our fit when perturbative evolution 
is assumed. The data points are values of extracted 
from individual simulation results for _R„ after extrapolating 
to a = and dividing out z(3, m^^) (n>4). Results are given 
for moments n = 4-10 and all 5 lattice spacings. Several points 
from different lattice spacings overlap in these plots. 



mass ratio: 

mb/rric = 4.53(4) (36) 

The other important output from our fit is a value for 
the parameter 

ao = aMs(5 GeV, = 3) = 0.2034(21). (37) 

To compare with other determinations of the coupling, 
we add vacuum polarization corrections from the c and 
b quarks, using the masses above, and evolve to the Z- 
meson mass p^| - |2l| : 

aj^{Mz,nf = 5) = 0.1183(7). (38) 



TABLE IV: Sources of uncertainty for the QCD coupling and 
mass determinations in this paper. In each case the uncer- 
tainty is given as a percentage of the final value. 







mi,(10) 


mb/rric 


mc(3) 


extrapolation 


0.2% 


0.6% 


0.5% 


0.2% 


perturbation theory 


0.5 


0.1 


0.5 


0.4 


statistical errors 


0.1 


0.3 


0.3 


0.2 


ruh extrapolation 


0.1 


0.1 


0.2 


0.0 


errors in ri 


0.2 


0.1 


0.1 


0.1 


errors in ri/a 


0.1 


0.3 


0.2 


0.1 


errors in m,,^ ,1^,^^ 


0.2 


0.1 


0.2 


0.0 


ceo prior 


0.1 


0.1 


0.1 


0.1 


gluon condensate 


0.0 


0.0 


0.0 


0.2 


Total 


0.6% 


0.7% 


0.8% 


0.6% 



Figure [2] shows how consistent our simulation results are 
with the theoretical curve for aj^{^, nf — 3) correspond- 
ing to our value for ag. For this figure we extracted 
values for from each i?„ separately by dividing out 
the dependence and z(3,to^^) using our best- fit pa- 
rameters, and then solving for ajj^ by matching with 
perturbation theory for r„. (In our fit, of course, we fit 
all i?„s simultaneously to obtain a single for all of 
them.) 

The dominant sources of error for our results are listed 
in Table IIVI The largest uncertainties come from: ex- 
trapolations to a — 0, especially for quantities involving 
b quarks; unknown higher-order terms in perturbation 
theory, especially for quantities involving c quarks; sta- 
tistical fiuctuations; extrapolations in the heavy quark 
mass, especially for quantities involving b quarks; and 
uncertainties in static-quark parameters ri/a and ri.The 
pattern of errors is as expected in each case. The non- 
perturbative contribution from the gluon condensate is 
negligible except for rric, again as expected; and errors 
due to mistuned sea-quark masses, finite volume errors, 
and uncertainties in MS coupling and mass evolution are 
negligible (<0.05%). 

The extrapolations of our data are not large. This 
is illustrated for rrih w nic in Figure [31 which shows the 
dependence of the reduced moments. The smallest 
two lattice spacings are sufficiently close to a = that the 
extrapolation is almost linear from those points. The a = 
extrapolated values we obtain here for the i?„ agree to 
within (smaller) errors with those in our previous paper: 
here we get 1.282(4), 1.527(4), 1.373(3), 1.304(2) with 
n — 4,6,8,10, respectively, for the masses used in the 
figure. 

We tested the stability of our analysis in several ways: 

• Vary perturbation theory: We chose fi = 3mh in 
order to keep scales large and (Xy^gdi) small. Our 
results are quite insensitive to fi, however. Choos- 
ing fi — mfi, for example, shifts none of our results 
by more than 0.2cr, and leaves all errors unchanged 
except for mc(3), where the error increases by a 
third. Taking /i = 9to?j shifts results by less than 
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FIG. 3: Lattice spacing dependence of Rn for masses m^^ 
within 5% of m,,^ and moments n = 4, 6, 8, 10. The dashed 
hnes show our fit for the average of these masses, and the 
points at a = are the continuum extrapolations of our data. 



0.4(7, and reduces the rric error by a third, leaving 
others only slightly reduced. Adding more terms 
to the perturbative expansions (Afpth = 6 — >■ 8) also 
has essentially no effect on the results. The prior 
for the unknown perturbative coefficients (Eq. (HI])) 
is twice as wide as suggested by our simulation re- 
sults (using the empirical Bayes criterion [13); we 
choose the larger width to be conservative. 

• Include more/fewer finite-a corrections: We set 
Nam — 30 for our results above. Using Nam — 15 
gives results that differ by less than 0.5a for nii, 
and much less for the other quantities. Much larger 
NamSi can be tested easily using the trick described 
in Section fill B 21 For example, replacing by 
i?^^" (Eq. (UHl)) with Nam = 80 and Nam = 30 
gives results that are essentially identical to those 
above. As discussed above, taking Nam = with 
the same Nam also gives the same results and is 
22 times faster (see the appendix for further dis- 
cussion). 

• Change n dependence of finite-a corrections: Re- 
placing the n-dependent prior for the expansion co- 
efficients (Eq. ([17])) by the n-independent prior 0± 
0.5 causes changes that are less than 0.3tT. The 
width of the original prior is optimal according to 
the empirical Bayes criterion — that is, it is the 
width suggested by the size of finite-a deviations 
observed in our simulation data. 

• Add more/fewer K/mrjh terms in z: Increasing the 
number of terms in the expansion for z from N^ — 4 
to 6 changes nothing by more than O.lcr. Decreas- 
ing to A^^ = 3 also has no effect. Again the width 
of the prior is optimal according to the empirical 
Bayes criterion. 

• Include more/fewer moments: Keeping all mo- 



ments 4 < n < 18 changes nothing by more 
than 0.5(7 and reduces errors slightly for every- 
thing other than m?,, where the errors are cut al- 
most in half: TOh(lO) = 3.623(15) GeV or mi,{mb) = 
4.170(13) GeV, both for 7i/ = 5. We continue to 
restrict ourselves to moments with n < 10 because 
these are the only moments for which we have ex- 
act third-order perturbation theory. Keeping just 
n = 4, 6 gives almost identical results for rric and 
q;^, with almost the same errors, but doubles the 
error on m^. 

• Omit simulation data: The coarsest two lattice 
spacings (configuration sets 1-5) affect our results 
only weakly. Leaving these out shifts no result 
by more than 0.5cr and leaves errors almost un- 
changed. Leaving out the smallest lattice spacing, 
however, increases errors significantly (almost dou- 
ble for Q!i;jg), while still shifting central values by 
less than 0.5(t. 

• Add large masses: Including cases with am^^ > 
1.95 from Table HI] leads to poor fits. The excluded 
data, however, do not deviate far from the best-fit 
lines. For example, the points marked with an "x" 
in Figure [1] are for the largest mass we studied, 
corresponding to m^^ =9. 15 GeV (last line in Ta- 

m^^ is too large for this case to 



ble|n|. Although 
be included in our fit, the values of i?n/r„ are only 
slightly below the fit results. 



V. NONPERTURBATIVE mt/m^ 

It is possible to extract the ratio of quark masses 
rrib/mc directly, without using the moments and with- 
out using perturbation theory. This provides an excellent 
nonperturbative check on our results from the moments. 

Ratios of quark masses are UV-cutoff independent and 
therefore the ratio of MS masses 



mi,{p,nf) _ mob 



rndp, Uf) 



moc 



(39) 



for any /i and nf, where mof, and toqc are the bare 
quark masses in the lattice quark action that give correct 
masses for the rjc and rjb, respectively. We obtain accu- 
rate mass ratios from this relationship by extrapolating 
to a — 0. We used such a method recently to determine 
mc/nis [Ul]. 

Here we have to modify our earlier method slightly 
because we cannot reach the 6-quark mass directly, but 
rather must simultaneously extrapolate to the b mass and 
the continuum limit. This is most simply done by deter- 
mining the functional dependence of the ratio 



w{m,j 



2moh 



m 



(40) 
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on the r]h mass and the lattice spacing. The ratio of 
MS masses is then given by the experimental masses of 
the T]c and T]b and the equation: 



(41) 



It might seem simpler to fit mp/i directly, rather than 
the ratio w] but using w significantly reduces the m^^ 
dependence (and therefore our extrapolation errors), and 
also makes our results quite insensitive to uncertainties 
in our values for the lattice spacing. 

We parameterize function w with an expansion mod- 
eled after the one we used to fit the moments: 



w{mrj^,a) ^ Zm{a) 1 + 




/ (42) 



Nam 



where, as for the moments, 



i+ i < max(A^t„„, N^). 



(43) 



Coefficients and Wn are determined by fitting function 
w{mjf^,a) to the values of 2amoh/(aw^^) from Table HIl 
The fit also determines the parameters Zm{a), one for 
each lattice spacing, which account for the running of 
the bare quark masses between different lattice spacings. 

The finite-a dependence is smaller here than for the 
moments, because the rjh is nonrelativistic (finite-a errors 
are suppressed by additional powers of v/c [1]), and the 
variation with m^^ stronger (twice that of z(3, rriri,^)). So 
here we use priors 



= ± 0.05 
w„ = ± 4 
Z„(a) = 1±0.5 



(44) 



with We again take iVam = 30, although identical 

results are obtained with Nam = 15. 

Our fit results are illustrated by Figure |4] which plots 
the ratio moh/nirj^ divided by moc/rrirj^ for a range of 
rjh masses. Our data for different lattice spacings is com- 
pared with our fit, and with the a = limit of our fit 
(sohd hue). The fit is excellent, with xV22 = 0.42 for 
the 22 pieces of data we fit (we again exclude cases with 
irj^ > 1.95). Using the ijc and rji, masses from Sec- 



tion |TVBl and Eq. (|1T|) with the best-fit values for the 
parameters, we obtain finally 



moc 



■ 4.49(4) 



as a- 



(45) 



which agrees well with our result from the moments 
(Eq. dSS])). 
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FIG. 4: Ratio moh/mr,^ divided by moc/m^^ (which we ap- 
proximate by w{m^^ 1 1)/2 from our fit) as a function of m,,^ . 
The sohd hne shows the a = extrapolation obtained from our 
fit. This is compared with simulation results for our 4 small- 
est lattice spacings, together with the best fits (dashed lines) 
corresponding to those lattice spacings. The point marked by 
an "x" is for the largest mass we tried (last line in Table [lT|) : 
this was not included in the fit because am,,^ is too large. 



VI. 



FROM WILSON LOOPS 



In a recent paper [20|, we presented a very accurate 
determination of the QCD coupling from simulation re- 
sults for Wilson loops. Here we want to compare those 
results to the value we obtain from heavy-quark corre- 
lators. First, however, we must update our earlier anal- 
ysis to take account of the new value for ri [13] given 
in Eq. (jlOp and improved values for ri /a fl3j given in Ta- 
ble m (The Wilson- loop paper uses some additional con- 
figuration sets: from Table II in that paper, sets 1, 6, 9, 
and 11 whose new ri/as are 1.813(8), 2.644(3), 5.281(8) 
and 5.283(8), respectively.) We have rerun our earlier 
analysis, updating ri, ri/a, and the c and b masses. The 
results are shown in Figure [SJ Combining results as in the 
earlier paper we obtain a final value from the Wilson-loop 
quantities of 



{Mz,nf = 5) = 0.1184(6), 



(46) 



with x^/22 = 0.3 for the 22 quantities in the figure. 
This agrees very well with the result in the earlier pa- 
per, ajj^{Mz) = 0.1183(8), but has a slightly smaller 
error, as expected given the smaller error in ri. This 
new value also agrees well with our very different de- 
termination from heavy-quark correlators (Eq. ([55)) V A 
breakdown of the error into its different sources can be 
found in Table IV of [l^ (reduce the ri and ri/a errors 
in that table by half to account for the improved values 
used here). 
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FIG. 5: Updated values for the 5-flavor at the Z-meson 
mass from each of 22 different short-distance quantities built 
from Wilson loops. The gray band indicates a composite av- 
erage, 0.1184(6). per data point is 0.3. 




FIG. 6: z{jj./mh,mr, 



(in GeV) for three different 



versus 

values of n/mh- The curve for /i = 3mh comes from the best 
fit to the moments. The other curves are obtained by evolving 
perturbatively from ^ = 3mh. 
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FIG. 7: Simulation results for reduced moments R„ with n = 
6, 8, 10 as functions of m,,^ for 5 different lattice spacings. 
The dashed lines show the corresponding behavior of our fit 
function, with the best-fit parameters. The curves for smaller 
lattice spacings extend further to the right. The solid lines 
show the a — limit of our best fit. 



VII. CONCLUSIONS 



In this paper, we improve significantly on our previous 
determinations of the QCD coupling and c-quark mass 
from heavy-quark correlators. This is principally due to 
the inclusion of a new, smaller lattice spacing in our anal- 
ysis. We also generated results for a variety of quark 
masses near rric, allowing us to interpolate more accu- 
rately to the physical value of mc- New third-order per- 
turbation theory makes i?io as useful now as i?4, Rq, and 
i?8 were in the earlier paper. Finally, in this paper, we 
fit multiple moments simultaneously, determining con- 
sistent values simultaneously for both the QCD coupling 
and the quark masses for all moments. Previously we ex- 
amined each moment or ratio of moments independently, 
extracting rricS or ajjgS independently of each other. Our 



new results, 

mc(3 GeV, n/ = 4) = 0.986(6) GeV (47) 
a^jg(Mz,n/ = 5) = 0.1183(7), 

agree well with our older results of 0.986(10) GeV and 
0.1174(12), respectively [l|. 

The much heavier b quark is usually analyzed using ef- 
fective field theories like NRQCD or the static-quark ap- 
proximation. By using very small lattice spacings and the 
very highly improved HISQ discretization for the heavy 
quarks, we are able to extend our analysis almost to the 
6-quark mass, using the same relativistic discretization 
that we use for c and lighter quarks. A 1.5% extrapo- 
lation of z{3,mh), from the largest m^^ used in our fits 
to rriri^, gives us a new, accurate determination of the 
&-quark mass, 

mb(10GeV,n/ = 5) = 3.617(25) GeV. (48) 
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This calculation demonstrates the utility of the HISQ 
formalism for studying b quarks on lattices that are com- 
putationally accessible today. This represents a break- 
through for b physics on the lattice since far greater pre- 
cision becomes possible when all quarks are treated us- 
ing the same formalism, and that formalism is relativis- 
tic and has a chiral symmetry. Even better would be 
to work right at the b mass, as opposed to extrapolat- 
ing from nearby; this would require a lattice spacing of 
order 0.03 fm. 

Both of our new c and b masses agree well with 
non-lattice determinations from vector-current correla- 
tors and experimental e+e~ collisions. A recent analysis 
of the continuum data gives 



mc(3 GeV, n/ 



4) 
5) 



0.986(13) GeV 
4.163(16) GeV 



(49) 



which compare well with our values of 0.986(6) GeV and 
4.164(23) GeV, respectively. This provides strong evi- 
dence that the different systematic errors in each calcu- 
lation are understood. 

Function z(/i/m/i, m,,^) is a by-product of our anal- 
ysis. It relates the MS quark mass m;j(/i) to the 
i]h mass (Eq. d?])). We show our result again in Figure [6] 
for ^ = 3mh, as well as for ^ = mh and ^ — which 
we obtain by evolving perturbatively from ii — Sruh- The 
latter two curves are relatively flat, and the last surpris- 
ingly close to 1 for most masses. 

Questions have been raised about the way perturba- 
tion theory is used in analyzing the perturbative parts 
of the moments [131 . Like [J we favor using larger 
scales than rric for c-quark correlators, but, as we have 
shown, our results are quite insensitive to /i over a broad 
range. Furthermore, the fact that our results, from 
pseudoscalar-density correlators, agree so well with the 
continuum results, from vector-current correlators, is also 
compelling evidence that perturbation theory is being 
handled correctly. We also find consistent results from 
several different moments, which is only possible if per- 
turbation theory is working well. Compare, for example. 
Figure [7] for the moments, as a function of m^^ , with 
the plots of Rn/vn in Figure [TJ Figure [7] shows very dif- 
ferent m,,^ behavior, at the 10-20% level, for different 
moments i?„; Figure [TJ where the perturbative part r„ 
is divided out, shows behavior that is almost moment- 
independent. 

An additional check on our use of perturbation theory 
comes from the close agreement between our perturbative 
result for the ratio mb/rric of MS masses (Eq. (|36p ) and 
our nonperturbative result for the ratio of HISQ masses 
(Eq. pS]) '). These should be and are equal to within 
our 1% errors. Taken together they suggest a composite 
result of: 



4.51(4) (composite). 



(50) 



Y decays 

X decays 

DIS [Fj] 

DIS [e,p ->jets] 

e"''e"[jets shps] 

electroweak 

e"''e"[jets shps] 

HPQCD: lattice wloops 

HPQCD: lattice current 
correlators 

World average: 
Bethke 0908.1 135 



0.11 0.115 0.12 0.125 
«s(Mz) 



0.13 



FIG. 8: The 5-flavor QCD coupling at the Z mass as 
determined by a variety of different methods. The non-lattice 
numbers used here are from the review in [2811. 



QCD coupling we get from the heavy-quark correlators, 
a^{Mz) = 0.1183(7), and that obtained from Wilson 
loops, 0.1184(6). These are radically different methods 
for determining the coupling. The first relies upon a con- 
tinuum quantity, extrapolated to a = 0, and continuum 
perturbation theory. The second relies upon quantities 
that are highly sensitive to the UV cutoff (vr/a) but are 
analyzed to all orders in the cutoff using lattice pertur- 
bation theory. Systematic errors are almost completely 
different in the two cases. The fact that they agree to 
within our 0.6% uncertainties is highly nontrivial evi- 
dence that perturbative and other potential errors are 
understood. 

Our coupling values also agree well with determina- 
tions from non-lattice methods. Figure [5] summarises re- 
cent results that were included in a world average by 
Bethke The world average resuff, 0.1184(7), was 

dominated by our previous determination from the Wil- 
son loop analysis. The average excluding our result was 
0.1186(11), which also agrees well. Including our new 
results into a new error-weighted world average gives 



(Mz) = 0.1184(4). 



The validity of our perturbative analyses is fur- 
ther supported by the close agreement between the 



Our new c mass is the most accurate currently avail- 
able. With it we can improve slightly on our recent de- 
termination of light quark masses using an accurate value 
for rric/ms, 11.85(16), derived completely nonperturba- 
tively from lattice calculations llj. Our new c mass, 
which becomes 1.093(6) GeV when converted to Uf — 2> 
at 2 GeV, implies: 

m^(2 GeV, n/ = 3) = 92.2(1.3) MeV, (51) 
md{2 GeV, = 3) = 4.77(15) MeV, 
m„(2GeV,n/ = 3) = 2.01(10) MeV. 

Our results for all 5 quark masses are compared with the 
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FIG. 9: MS masses, for the 5 lightest quarks, from this pa- 
per compared with the Particle Data Group's current esti- 
mates [231 . Each mass is quoted at its conventional scale: 
2 GeV for it, d, s (n/ =3); rric for c (n/ =4); irib for b (n/ = 5). 



Particle Data Group's 2009 values in Figure [9l Agree- 
ment is excellent, but our uncertainties are much smaller 
in every case, and by an order of magnitude for the 
strange and light quarks. 

Finally we note that the consistency between quark 
masses from lattice and non-lattice analyses, and be- 
tween couplings from heavy-quark correlators and Wil- 
son loops provides further evidence that taste-changing 
interactions in the HISQ and ASQTAD quark formalisms 
are understood and vanish as a -> 0. While early con- 
cerns about the validity of these formalisms have been 
largely addressed both by formal arguments Il3l. I30i434| 
and by extensive empirical studies [8l-[lll [26l. l35[ - l39l |. it 
remains important to test the simulation technology of 
lattice QCD with increasing precision, given the growing 
importance of lattice results for phenomenology. 
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In Section llll B 2l we used a trick to simplify our fits by, 
in effect, transferring fit terms from the fit function into 
the errors of the fit data. This trick can greatly speed 
up complicated fits. Here we present a formal deriva- 
tion of this procedure for three increasingly complicated 
situations. 



A. Linear Least Squares — Exact Data 

Assuming we know D values jji for a quantity y which 
can be expressed as a power series in x, 



E 



(52) 



we wish to obtain a best fit for the first F unknown coef- 
ficients c„. The c„ are then our random variables. If we 
are able to make reasonable estimates for their means and 
standard deviations (t„, in the absence of additional in- 
formation, maximizing entropy suggests a Gaussian prior 
of 



P{c) oc e 



(53) 



For simplicity, we assume throughout that the c„ are 
uncorrelated and have a prior mean of zero; extending to 
more general cases is straightforward. 

If we knew all coefficient values, then the data yi would 
be completely determined, with 



D-l 



p{y\c)^ Yl %,-^c„xn 

i—0 n 

Bayes' theorem 

P{c\y) cx P{y\c)P{c) 



(54) 



(55) 



allows us to convert this into a distribution for c given 
the data y. 

If we are only interested in fitting a subset of coeffi- 
cients with n < F, we integrate over the remaining 
Cn> , giving 



P(c<|?/)oce "X 



(56) 



c„a;" e " 



We replace the delta function by its Fourier representa- 
tion, integrate over first the , then the Fourier vari- 
ables, to obtain 



P(c<|y)cxe "X 
(detai)- 



(57) 
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Here 

Ay^ EE - ^Cn^x" (58) 

is the discrepancy between the measured yi and the por- 
tion of the series to be kept in the fit, the dot product 
sums over the D data points, and 

n> 

The correlation matrix <t\ is independent of c< (so the 
determinant is constant), and is the same as one would 
compute directly by 

(Ay, Aj/,)^^ = r£ xT J2 ^" ) (60) 

using 

(Cm> C„>)^_^ = CT,^, Smn ■ (61) 

Finally, we fit c< by minimizing x^, which includes these 
correlations and is augmented by the remaining c< pri- 
ors. Because the distribution is Gaussian, the c< at their 
minima are equal to their average values. 

The correlation matrix a\ properly accounts for cor- 
relations in the discrepancy, due to the neglected terms, 
between y and the portion of the series retained. If F 
terms are kept in the series, cta is 0{x^), enforcing agree- 
ment between y and the finite series to this order, as 
appropriate. It also suggests an alternative but equiva- 
lent approach. We may define new random (rather than 
exact) versions of y, whose correlation matrix is ct^, by 
moving the c> terms to the left side of Eq. ([5^ . Using 
the truncated series as a model for these random data, 
straightforward application of Bayes' theorem [13] again 
implies the distribution in Eq. (j57p . 

One useful consequence is that, as long as we include 
the correlations for the c> , we may arbitrarily reduce the 
number of coefficients c< retained, even to as few as one, 
and still obtain the same minimization values. To see 
this, note that to compute a particular (c„^), we could 
start with the full distribution and integrate over all cs. 
The integral over c> produces P(c<|?/), which we then 
use in the integral over c<; the result will be the same 
regardless of where the dividing line is set, as long as 
it does not include c„^ . (We could even include in ct^ 
terms of order less than n.) Because averaging and min- 
imization give the same result, the minimization value 
for c„^ will also remain unchanged. This is also true of 
the c„^ error. While the result is the same, reducing 
the number of terms in the series to fit can significantly 
improve the fitting time. 

B. Fits to Nonlinear Functions — Exact Data 

We now consider fitting to the data yi a general func- 
tion gi{cn) not necessarily linear in the parameters c, and 



where we assume yi = 5i(c„) exactly for properly chosen 
c„. Now 

D-l 

P{y\c)<x\{5{y,-Y,g,{cn)). (62) 

■i— n 

Combining with the prior P{c) and integrating over the 
c> gives P(c<|y). 

If our estimate of prior means is good, expanding g 
around c> = should give a reasonable approximation; 
an expansion to first order gives a Gaussian. More specif- 
ically, defining 

ffi(c<) = 5j(c<,c> = 0) (63) 

and 

Ay, = -5,(c<) , (64) 

and integrating over c> in this Gaussian approximation 
gives as before 

P(c< |y) oc e"^"< ''"</^'""x (65) 
(detai(c<))-i/^ g-A,.(2.i(c<))-i.Aj,^ 

but with 

o-A (c< )»j = X! ^"5' ('^< ) ^« ("^< ) • (66) 

This is again the correlation one would compute directly 
for (Ayi Ayj)^^ after expanding g to first order in c>. 

We have not expanded in c< , so (j\ depends on c< , the 
determinant in front is not constant, and the dependence 
of Ay on c< is not in general linear. In practice, how- 
ever, we will often further approximate the distribution 
by setting the c< to their prior means in cr^(c<) before 
minimization. 

Because 5(c<) is nonlinear, c< from minimization can 
differ slightly from (c<), and due to approximations 
made, can vary somewhat with the number of terms re- 
tained. 



C. Fits to Data with Intrinsic Statistical Errors 

Finally we consider the most general case, in which 
the data y contribute intrinsic statistical uncertainties in 
addition to those associated with the truncated series. 
If we measure a range of values for y with an average 
(y) and correlation matrix (7^, then for sufficiently large 
samples we expect a Gaussian distribution 

P((y) |c) oc e-(<'^>-f('=»-(2'^«)"'-(<«>-ff('=» (67) 

rather than the delta function above. Combining with 
the prior P(c) gives P(c| (y)). 
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Expanding .9i(c<,c>) to first order around c> = 0, 
defining 

Ay, = (y,) -5.(c<) , (68) 
and integrating P{c\ (y)) over c> gives 

P(c<l(y)) oce"^"< X (69) 

(det a2^(c<))-i/2 g-A,.(2.2,(c<))-.A,^ 

The resulting correlation matrix is a combination of true 
statistical and neglected series contributions, with 

al^ic<) = al + al{c<) , (70) 

as one would obtain by including both sources of uncer- 
tainty in computing (Ay^ Ay^) directly. With no sta- 
tistical fluctuations in y, cr^ = 0, and it reduces to the 
previous result. When cr^ is nonzero but small, a\ still 
makes an important contribution. 



D. Application to this Paper 



We used the technique described here in much of our 
testing and tuning (but not for our final results) to speed 
up the (am^^/2)^ fit. As described in Section [ill B 21 
we kept corrections through order Nam — 80 but moved 
all but Nam ^ Nam out of thc fit fuuctiou and into the 
errors for the reduced moments. If we set Nam — 3, for 
example, our fit to the i?„ simulation data changes from 
Figure [7] to Figure ITOl The small Nam means that each 
point in Figure [TU] has much larger error bars, coming 
from terms moved into the i?„s. The final fit 

results, however, are almost identical in both cases (to 
within less than O.Ict), with the same errors. Note that 
the Rn errors in Figure [TU] are highly correlated, which is 
why the fit curve passes through the central value for each 
point. As discussed above these correlations are essential 
if results are to be independent of the value of Nam ■ 
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